function [part] = partial_LogGMM(beta2, bpartial, mean, mean_l, mean_h, U, RAND, n, k)
global draws

% This file estimates the partial effects at the mean (zero for binary variables)
% Tobias Pfutze July 2009


% Compute matrix of drawn errors ER
for j=1:draws
    PART(1:k-1,j)=(bpartial(1:k-1)*exp(mean*beta2(1:k)+U(RAND(j),1))...
        +(bpartial(1:k-1)-bpartial(k:2*k-2))*exp(mean*(beta2(1:k)+beta2(k+1:2*k))+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean*beta2(1:k)+U(RAND(j),1))+exp(mean*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    PART(k:2*k-2,j)=(bpartial(k:2*k-2)*exp(mean*beta2(k+1:2*k)+U(RAND(j),2))...
        +(bpartial(k:2*k-2)-bpartial(1:k-1))*exp(mean*(beta2(1:k)+beta2(k+1:2*k))+U(RAND(j),1)+U(RAND(j),2)))...
        /((1+exp(mean*beta2(1:k)+U(RAND(j),1))+exp(mean*beta2(k+1:2*k)+U(RAND(j),2)))^2);
    
    %Replace discrete changes for binary variables
    
    PART(3:5,j)=exp(mean_h*beta2(1:k)+U(RAND(j),1))...
        ./(1+exp(mean_h*beta2(1:k)+U(RAND(j),1))+exp(mean_h*beta2(k+1:2*k)+U(RAND(j),2)))...
        -exp(mean_l*beta2(1:k)+U(RAND(j),1))...
        ./(1+exp(mean_l*beta2(1:k)+U(RAND(j),1))+exp(mean_l*beta2(k+1:2*k)+U(RAND(j),2)));
    
    PART(k+2:k+4,j)=exp(mean_h*beta2(k+1:2*k)+U(RAND(j),2))...
        ./(1+exp(mean_h*beta2(1:k)+U(RAND(j),1))+exp(mean_h*beta2(k+1:2*k)+U(RAND(j),2)))...
        -exp(mean_l*beta2(k+1:2*k)+U(RAND(j),2))...
        ./(1+exp(mean_l*beta2(1:k)+U(RAND(j),1))+exp(mean_l*beta2(k+1:2*k)+U(RAND(j),2)));
end

part=PART*ones(draws,1)/draws;
